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A simple solid-on-solid model of adsorbate - induced faceting is studied by using a modified 
Wang-Landau method. The phase diagram for this system is constructed by computing the density 
of states in a special two-dimensional energy space. A finite-size scaling analysis of transition 
temperature and specific heat shows that faceting transition is the first order phase transition. 
Logarithmic dependence of the mean-square width of the surface on system size indicates that 
surface is rough above the transition temperature. 



PACS numbers: 68.35.Rh, 68.43.De, 64.60.Cn 



I. INTRODUCTION 

It has been recently demonstrated [iHl] that surfaces 
such as W(lll) and Mo(lll) covered with a single phys- 
ical monolayer of certain metals (Pd, Rh, Ir, Pt, Au) 
undergo massive reconstruction from a planar morphol- 
ogy to a faceted surface upon annealing at T > TOOK. 
The faceted surface is covered by three-sided pyramids 
with mainly {211} facets. A major mechanism for facets 
formation is minimization of total surface free energy [6]. 
Faceting of bcc(lll) and fcc(210) surfaces can also be 
induced by oxygen or other nonmetallic impurities (see 
Ref. 7] and references therein). Adsorbate - induced 
faceting is also observed on curved surfaces and 
this phenomenon has been recently used in fabrication of 
electron and ion point sources [Hi [l3 • 

Investigations of thermal stability of faceted sur- 
faces have revealed that reversible phase transition 
occurs in systems Pd/Mo(lll)[i|, 0/Mo(lll)[i|, and 
Pd/W(lll)[a|- As a result of faceting transition, the 
faceted surface changes into planar surface at the transi- 
tion temperature. 

In theoretical studies of complicated surface problems 
(e.g. roughening transition, surface reconstruction, sur- 
face growth, surface phase transitions) simple models 
like lattice gas models or solid-on-solid (SOS) models 
are applied |l"3l-[l7| . In our earlier paper JlS*] we have 
proposed a SOS model to study the adsorbate-induced 
faceting of the bcc(lll) crystal surface at constant cov- 
erage. Monte Carlo simulation results show formation 
of pyramidal facets in accordance with experimental ob- 
servations. Moreover, the model describes a reversible 
phase transition from a faceted surface to a disordered 
(111) surface. This model reproduces also formation of 
{211} step-like facets on curved surfaces [Toj. 



Simulation results [18| - based on the Metropolis algo- 



rithm, indicate that the faceting is the first-order tran- 
sition. However, identification of the nature of this tran- 
sition in finite-sized system is difficult and it can be 
solved by use of finite-size scaling [23|. Another way to 
identify the first-order transition is to study the distri- 
bution of energy which has a double-peak structure in 
the vicinity of transition temperature. This can be eas- 
ily accomplished by use of the Wang - Landau (WL) 



method|2ll |22[. The WL method is based on accurate 
calculation of density of states and therefore allows for 
calculating thermodynamical functions. The method is 
esp ecially useful in investigation of phase transitions pTI - 
|24| to determine the order of transition, transition tem- 
perature and behavior of thermodynamical quantities. 

In this paper we use WL method to study faceting 
transition in over layer-induced faceting on bcc(lll) sur- 
face. A short presentation of SOS model of adsorbate - 
induced faceting is given in Sec. |lTl Application of WL 
method and its modification proposed by Belardinelli and 
Pereyra[25] to calculation density of states is presented 
in Sec. IIIIl Sec. IIVI contains a method of construction 
of a phase diagram for a finite system with competitive 
interactions. In this method, energy of the system is 
decomposed into a few parts and density of states is cal- 
culated in multi-dimensional energy space. Due to this 
method, the phase diagram is constructed for the SOS 
model of finite sizes. In order to determine the order 
of phase transition, the finite - size scaling of the transi- 
tion temperature and specific heat is presented in SecfVl 
The nature of high-temperature phase is investigated in 
Sec. IVIl Contrary to experimental results, the finite size 
analysis of the mean-square width of the surface indi- 
cates that above the transition temperature a surface is 
rough. 



II. THE SOS MODEL 

To study an adsorbate-induced faceting on bcc(lll) 
surface a solid-on-solid model has been proposed [l8|. 
The model consist of columns placed on the triangular 
lattice obtained by projection of the bcc crystal lattice 
on the (111) plane (see Fig. [T]). A column height hi at 
site i = (ixjiy) in the Ith sublattice, I = 0,1,2 takes 
discrete values of the form hi — Sui + I where Ui is the 
number of atoms in this column. The model is designed 
for constant coverage of 1 physical monolayer - critical 
coverage for adsorbate-induced faceting This means 
that each column has exactly one adsorbate atom placed 
at the highest position. There is also a restriction im- 
posed on column heights: the nearest neighbor column 
heights can differ only by ±1, ±2. 
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FIG. 1. Schematic top view of the bcc(lll) surface. Atoms 
from three successive geometrical layers represent positions of 
columns in the SOS model. The Z axis is normal to the (111) 
plane and ax — a^/2/2, ay = a\/Q/i. 

The surface formation energy can be expressed as in- 
teraction energy of columns 

« [ ii 

+ E [2^20..,,, (0) + (2 J2 + K2) c^..,, (3)] (1) 

+ J2E N J3 (2) + ^« j3 (4)] \+NJo, 

where the sums over ji , j2 , and J3 represent the sums over 
first, second, and third neighbours of the column at site 
«, respectively. The function ujij{k) — 5 {\hi — — fc), 
for k — 0,1,2,..., is expressed by the Kronecker delta: 
5{x) = 1 for a; = and 5{x) = for x ^ 0. Model pa- 
rameters Jo, Ji, J2, Ki, and K2 depend on interaction 
energies between substrate and adsorbate atoms (for de- 
tails see [III)- 

It turns out that the energy of the nearest neighbours 
interaction in the Hamiltonian is conserved. This prop- 
erty follows from symmetry of the model and the restric- 
tion imposed on the on column heights. Thus the en- 
ergy of the nearest neighbours interactions can be treated 
as the reference energy. Choosing coupling constant J2 
as the unit of energy, we will work with one model pa- 
rameter K ~ K2I Ji- The dimensionless Hamiltonian 
% = %/ Ji, takes the following form 

^ = lE|E[2w.,..(0) + (2 + i^)u;.,,,(3)] 

i I 72 

^ . (2) 

+ EKj3(2)+^.j3(44+'^,ef, 



In what follows we will use dimensionless energy omitting 
the tilde and "Href- It has been shown [IQ] that the energy 
of the (211) face is minimal when — 2 < if < 0, whereas 
the (111) surface is stable for i<r > 0. For K < —2 the 
(Oil) face has minimal energy. 



III. DENSITY OF STATES 

In this paper we calculate density of state for the SOS 
model using a modification of the Wang-Landau method 
[2^ 1 proposed by Belardinelli and Perevrapsl l26l|. 



A. The Wang — Landau method 

The WL method is based on a random walk which 
produces a fiat histogram in the energy space. A trial 
configuration of energy E' is accepted with probability 

W{E ^ E') ^min[l,g{E)/g{E')], (3) 

where E is energy of the current configuration and g{E) 
means the density of states. However, the g{E) is not 
known and calculation of the density of states is the 
main goal of WL method. To achieve this goal, the 
g{E) is changed after each step of the random walk 
g{E) — )■ fg{E) by a modification factor / > 1. More- 
over, it is assumed that initially g{E) = 1 for all energies 
and / = e. 

Very recently, Belardinelli and Pereyra have demon- 
strated [2^ saturation of errors in the WL method, or 
nonconvergence of calculated density of states to the ex- 
act value. Moreover, they have shown that if the refine- 
ment parameter ln/(i) depends on time as oc for 
large time, the calculated density of states approaches 
asymptotically to the exact values as oc This fact 

is used in Belardinelli - Pereyra (BP) modification of the 
WL method. 

To present the BP method let us introduce some quan- 
tities. It is assumed that the random walk is performed in 
the energy range Emin ^ E l£ Emax with f2 different en- 
ergy levels. A Monte Carlo time t = j /Vl\s defined as the 
number of trial configurations j used so far with respect 
to the number of energy states 51. From numerical reason 
it is convenient to use S{E) = \\ig{E) and F = \\\f in- 
stead of g{E) and /. There are two stages in calculation 
of S{E) in the BP method. In the first stage, the refine- 
ment parameter F is changed similarly as in the original 
WL method, i.e., F F/2 when all energy states are 
visited in the random walk with given F. Please notice 
that the criterion for fiatness of the histogram is not used 
here. The second stage begins at critical time tc defined 
as a moment when the new value of F becomes smaller 
than 1 /t. From this time, the refinement parameter takes 
the continuous form F{t) — t^^. The second stage lasts 
until the refinement parameter reaches a predefined value 
-Fmin (typically Fmin = 10~*). To control the convergence 
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of S(E,t) during the random walk it is useful to measure 
the following quantities: the histogram H{E, t) and its 
averaged value at time t, {H{t)) = X^b ^(-^i ^^'^ 
the width of the histogram AH{t) — -ffmax(i) — ^min(^), 
and the relative width SH{t) = j^^- According to re- 
sult of Ref.jlBl the quantity SH{t) has the same longtime 
behaviour as the error of S{E, t) 



5H{t) oc r 



(4) 



As the exact value of S{E) is not known, the time de- 
pendence of SH(t) can be used to evaluate convergence 
and the accuracy of BP algorithm. 



B. Calculation of S(E) for the SOS model 



moments of energy 



J2 E"^ eKv{S{E) - E/T) 

_E 

j:eMS{E)-E/T) ' 

E 



(6) 



the specific heat per site 

{E^)-{E)^ 

the Binder's fourth cumulant 



(7) 



V4 = l- 



3{E^y 



(8) 



We consider the SOS model on the rectangular lat- 
tice with Nx and Ny columns along x and y axis, re- 
spectively, and with periodic boundary conditions. A 
relation = \Ny is assumed to assure approximate 
equality of linear lattice sizes along the x and y axis. 
The linear system size is defined as L = y'N^Ny. The 
quantity S{E) = hig{E) is calculated by performing the 
the random walk in the energy space with the tran- 
sition probability given by Eq. ([3]). At each step a 
trial configuration is generated by choosing two lattice 
sites i and j and changing heights of columns at these 
sites: {hi, hj) — >■ {hi — 3, hj + 3). Due to constrains im- 
posed on columns height in the SOS model, the change 
{hi,hj) — {hi — 3,hj + 3) is allowed only when hi is a 
local maximum and hj is a local minimum. Local max- 
imum (minimum) at site k denotes that column hk is 
higher (lower) than its 6 nearest neighbor columns, re- 
spectively. In order to speed up calculation we use lists 
of local maxima and minima and a trial configuration is 
generated by random choice of a maximum and a mini- 
mum. 

During the preliminary application of the BP method 
to the SOS model we encountered a problem of very large 
tc even for the small linear sizes of the system. Contrary 
to simulation of Ising model, the number of steps needed 
to visit each energy level at least once becomes very large 
even for the first value of the refinement parameter _Fo = 
1. To overcome this problem, we modify the first stage 
of the BP method by introducing a separation s between 
successive increments of S(E) and H(E) in the random 
walk, similarly as in Ref.[27]. 

The errors SS{E) of S{E) are estimated from sample 
of m independent Si{E), . . . Sm{E) results of simulations 

_ Y.^'S(E) 



and then an average error 5S 



is calculated. 



Having calculated S{E) one can easily investigate tem- 
perature dependence of various quantities discussed in 
the paper: 

the energy distribution 



P{E,T) 



cxp(5(i;) - E/T) 
Y.eMS{E)-E/TY 



(5) 



IV. PHASE DIAGRAM 

WL and BP methods can be easily applied to con- 
struct a phase diagram for a finite system with com- 
petitive interactions by performing the random walk in 
multi-dimensional energy space. In this paper we con- 
struct the phase diagram in {T,K) plane by use of a 
modified BP method and performing the random walk 
in two-dimensional energy space {Ej,Ex)- To do this 
the energy of the system, E = H — "Hrct from Eq. is 
decomposed into two parts E — Ej + KEk, where KEk 
represents interaction energy with the coupling constant 
K and Ej stands for remaining contribution to E. 

Having calculated density of states g{Ej,EK) or 
S{Ej,Ek) = \ng{Ej, Ek) one can easily obtain the 
mean energy (E), the specific heat C, and other quanti- 
ties for any value of the coupling constant K. For exam- 
ple, the mean energy can be calculated as 



{E{K)) = 



Y.Y.{E., + KEK)e^^'^-'-^'<y 

E.J Ek 



Ej Ek 



S{E.,,Ef, 



(9) 



This approach is limited to rather small systems because 
the number of states n,{Ej, Ek) is much more greater 
than fl{E) in one-dimensional energy space. For exam- 
ple, Q.{Ej,Ek) reaches 82073 and 464162 for L = 26 
and 39, respectively whereas Q,{E) amounts to 304 and 
718 for K — —1. Therefore, we limit study of the phase 
diagram to two system sizes, L = 26 and 39. 

In simulations of S{Ej, Ek) we used separation s — AL 
and the refinement parameter F was limited by -Fmin = 
10~^. For such parameters and system size L = 39 the 
computation takes about 15 days on a 2.6 GHz Opteron 
processor. On the other hand, the average error of 
S{Ej, Ek) estimated on results of four independent runs, 
was rather small: 3 * 10""^ and 2 * 10"'^ for L = 26 and 
L = 39, respectively. Hence, this validates the use of 
S{Ej, Ek) to calculate thermodynamical quantities. 
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FIG. 2. Phase diagram for system size L = 39. Labels D, 
F112, Foil, and M denotes disordered, faceted {112}, faceted 
{Oil}, and mixed phases, respectively. 



To construct the phase diagram in the (T, K) plane 
for a finite system of a linear size L we treat temperature 
Tc{L,K) at which the specific heat has a maximum, as 
temperature of a phase transition. A simple optimiza- 
tion method - golden section search [2^ was applied to 
locate the maximum of the specific heat as a function of 
temperature for a given K. We also examined tempera- 
ture dependence of the Binder's fourth cumulant Eq. 
because V4 has a minimum at a first-order phase transi- 
tion. 

The phase diagram (see Fig. [2]) constructed for the 
lattice size L — 39 comprises three phases: faceted 112, 
faceted Oil, and disordered. There is also a region be- 
tween two faceted phases where the mixture of these 
phases appears. The line separated a faceted phase and 
the mixed one is determined from location of the ad- 
ditional maximum in the specific heat (see for example 
Fig. [31). Thus, in this case the first peak corresponds to 
transition from faceted 112 phase to the mix one, whereas 
the second peak is generated by transition from faceted 
112 phase to the disordered phase. The phase diagram for 
a smaller linear size L = 26 has the same qualitative form 
as the diagram in Fig.[5]but the phase transition lines are 
slightly shifted due to the finite-size effects. The largest 
differences are observed for transition temperature from 
the F112 phase to the disordered one (see Fig. 2]). In 
what follows we limit our consideration to this transi- 
tion because it is observed experimentally [l|, @ . Big dif- 
ferences between transition temperatures for these small 
systems require the finite-size analysis for larger linear 
system sizes. This can be performed by calculation of 
density of states in one-dimensional energy space. In- 
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FIG. 3. Specific heat with two peaks for K = —1.90. 

spection of energy distribution indicates that transition 
between the F112 and the disordered phase is of first 
order phase transition because the energy distribution 
at transition temperature has the double peak form (see 
Fig. [5]). The first peak represents the F112 phase and the 
second peak represents the disordered phase. The nature 
of this phase transition is also confirmed by temperature 
dependence of the Binder's fourth cumulant V4, which 
has minimum at temperature close to Tc{L). Thus, it is 
enough to perform a detailed analysis for one value of the 
coupling constant K , and we choose the case K ^ —1 io 
minimize the number of energy states or computing time. 



V. FINITE-SIZE SCALING 

The S{E) is calculated for the coupling constant K = 
— 1 in a single energy interval (i?min, -E^max) in order to 
avoid boundaries errors caused by partition of energy 
space into several pieces. It is especially important in 
a case of first order phase transition, where at transition 
temperature Tc, the energy distribution Eq. ([S]) has two 
peaks of equal heights ai E = E- and E = E+. In our 
case the interval {E-,E+) makes up about 1/3 of the 
whole energy interval. The energy of planar face (111) 
is chosen as i?max and this state has known degeneracy 
g{Eraax) = 1- On the other hand, it is difficult to reach 
states close to the minimal energy due to edge energies, 
hence we choose -Emin as small as possible to assure con- 
vergence of S{E). 

To study the size dependence of some physical quanti- 
ties the following numbers L where used in calculations 
26, 39, 52, 65, 78, 91, and 104. For each value of L 
the S{E) was calculated until refinement parameter F 
reached the minimal value Fn^n = 10~®. Values of sepa- 
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FIG. 4. Finite-size effect. The transition line between the 
faceted F112 and disordered D phases. Continuous (dashed) 
line represent system size L 39 (26), respectively. 
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FIG. 5. Energy distribution at transition temperature for 
three values of K and system size L = 39. 



ration s used in the first stage of BP method were com- 
parable with the number of energy states (s w 0.7J7). To 
estimate errors of S{E) each calculation were repeated 
5 times and average density was obtained. The relative 
averaged error of S{E) was smaller than 3 x 10^'* for 
L < 104 and 1.4 x 10~^ for L = 104. The increase of 
the error of S{E) for the largest size studied here was 
caused by large fluctuation of the histogram in the low 
energy range. Therefore, we did not study the systems 



FIG. 6. A semi-log plot of energy distribution at transition 
temperature for 3 linear sizes of system. 



with sizes L > 104. 



A. Scaling of transition temperature 

We calculated Tc using fact, that the energy dis- 
tribution Eq. ( [5|) has two peaks of equal heights at 
T = Tc. The results presented in Fig. |6] demonstrate 
that P{E, Tc) has two maxima: the first at energy E^ 
of faceted phase and the second at energy E+ of dis- 
ordered phase. The minimum between these peaks be- 
comes deeper as system sizes increases. The probability 
to find the system at minimum P(£'rnin) is 7 x 10~^, 
2 X 10"^ and 3 x 10"^" of PiE±) for system size L 
52, 78, and 104 respectively. From energy distributions 
at T — Tc{L) one can calculate the free-energy barrier 
AT{L) = ~Tc{lnP{E^in,L)-\nP{E_,L)) which should 
scale as AT{L) ~ L"*"^ at a first order transition 
This scaling is confirmed in our calculation because for 
L > 50 the free-energy barrier has a linear form fitted 
by AJ"(L) = 0.354L - 13.01. 

Having calculated Tc{L) for several linear system sizes 
we study scaling of Tc{L) to calculate the transition tem- 
perature in the limit L — >■ oo. Our results agree with 
theory of scaling at first order phase transition 



TciL) = Tc 



aiL ^ 



(10) 



The second term proportional to L^* is needed for sys- 
tem size L < 75 (see Fig. [7]). The transition temperature 
Tc — 1.12852 is obtained from fitting the results with 
L > 50 by the function from Eq. (fTO)l . Similar result 
Tc ~ 1.12829 is obtained by applying the Burlisch -Stoer 
extrapolation [28| to Tc{L). We also studied scaling of 
temperature T^[L) at which the specific heat has a max- 
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FIG. 7. Transition temperature versus . Errors are 
smaller than the symbol size. Thin and thick lines denote 
linear (in L~^) and non-linear fitting of simulation result by 
the scaling function (Eq. HOjl). 



FIG. 8. Scaling of the specific heat. Line represents the 
scaling function (Eq. symbols represent simulation data 

and label denotes linear system size L. 



imum in system of linear size L. The differences between 
Tc(L) and T^{L) are of order 10~^ and extrapolation of 
T^{L) for L ^ oo yields = 1.12837. From these two 
results we estimate transition temperature for K = —1 
as Tr = 1.12844 ±0.00008. 



B. Scaling of specific heat 

We found that scaling of the specific heat of the sys- 
tem discussed in this paper can be well described by the 
function of the form used by Challa et al. [29] for q-state 
Potts model. 



C{L) = 



[exp{X{L)) + DeMX{L))] 



(11) 



where 



and 



Q{L) 



E+{L)-E^{L) 



X{L) 



2T 



The parameter D in Eq. (jlll) replaces the expression 
q(C_/C+)^/^ for the q-state Potts model. The value 
of 13 = 1.62 minimizes the sum of deviation of simulated 
data from the scahng function for L > 60. 

As seen in Fig. |S1 the scaling function (Eq. ([TT])) well 
describes the shape of the specific heat near transition 
temperature Tc{L) for L > 60. For smaller systems stud- 
ied here ( L < 60 ) this scaling does not apply - it yields 
incorrect position and height of specific heat maximum. 



VI. DISORDERED PHASE 



A nature of disordered phase is not clear. LEED exper- 
iment results [H suggested existence of flat phase above 
the faceting temperature. On the other hand, previous 
MC simulation[l^ indicated that this phase is not fiat 
one but disordered faceted phase characterized by chaotic 
hill-and-valley structure which comprises of randomly 
distributed small facets mainly of {112} orientation. To 
clarify this problem we study the size-dependence of the 
mean-square width of the surface 



6h^ 



(12) 



where sum runs over lattice sites and h is the arithmetic 
average of column heights. This quantity has been used 
in investigation of roughening transition [T^ - [T6| in SOS 
models because Sh^ as function of a system size L has 
logarithmic dependence in the rough phase. 

In order to investigate the dependence of Sh^ on tem- 
perature and system size L we computed the microcanon- 
ical averages Sh^ (E) for each energy state in systems with 
L < 105. The microcanonical averages Sh^{E) were cal- 
culated by performing the random walk in the energy 
space using earlier computed densities of states g{E) via 
BP method. To assure high accuracy of 6h^{E) each 
state E was visited on average 10^ times. 
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FIG. 9. Semi-log plot of 5h^ versus T for several lattice sizes 
L. 



The canonical average 5h?'{T) is obtained as 



J2Sh%E)exp{SiE) - E/T) 

_E 

j:exp{S{E)-E/T) 

E 



(13) 



This way of calculating 5h?{T) is similar to a method of 
computing moments of magnetization [30| . 

Simulation results show (see Fig. [5]) that Sh^ (T) is de- 
creasing function of temperature and it rapidly changes 
at transition temperature Tc{L). On the other hand, Sh^ 
is increasing function of linear system size L. In the high- 
temperature phase dh^ scales logarithmically with L in 
the whole temperature range T > Tc (see Fig. [TU)) and 
results of simulations are very well fitted by the function 



6h^{T,L) = A{T)\nL + B{T) 



(14) 



The amplitude A has the largest value at T = Tc and it 
decreases down to a saturation value 0.4625 for T > Tc 
(see Fig.[TT]). The logarithmic dependence of the mean- 
square width of the surface on linear size of the system 
indicates that the disordered phase is a rough phase. 
Hence, the faceting transition is the first-order rough- 
ening phase transition. On the other hand, the surface 
above Tc, has a disordered hill-and- valley structure as 
follows from Monte Carlo simulations (see e.g. Fig. 5 
in Ref.fisl'). Hence, this is different type of rough sur- 
face than that observed in typical roughening transition 
where a surface becomes rough by formation of steps [l6|. 



VII. DISCUSSION 

It is demonstrated that the Wang-Landau method can 
be applied to construct a phase diagram for a system 
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FIG. 10. Semi-log plot of size-dependence of Sh'^ at several 
temperatures T. 




FIG. 11. Plot of amplitude A versus temperature T. 



with competitive interactions if one performs calculations 
of density of states in a special two-dimensional energy 
space. From such density of states one can compute ther- 
modynamical quantities for any value of interaction ener- 
gies. Then one can find phase transition lines and nature 
of transitions. However such approach is limited to rather 
small linear sizes ( in our case for L < 40) due to huge 
numbers of states occurring in two-dimensional energy 
space. Because of that, a finite-size analysis should be 
performed for chosen values of interaction energies to find 
the order of the phase transition, transition temperature 
and other interesting quantities. 

We found that transition from faceted {112} phase to 
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high-temperature phase is of first order because the en- 
ergy distribution has double-peak form at transition tem- 
perature. The transition temperature scales as L^^ for 
large linear system size L, but higher order corrections 
are needed to the scaling for L < 75. Scaling of the spe- 
cific heat is well described by the function of the form 
used by Challa et al. ^2§\ for q-state Potts model. 

In order to clarify the nature of high-temperature 
phase we calculated the mean-square width of the sur- 
face Sh^ . Experimental results [l| suggest that surface 
becomes flat above transition temperature. However, in- 
vestigation of surface morphology in high temperatures 
is not so easy. In most of experiments, such investiga- 



tions are performed after quick cooling the sample to low 
temperatures. In case of adsorbate induced faceting, at- 
tempts of freezing of high-temperature phase were always 
failed for available cooling rates fl'|. On the other hand, 
the nature of high-temperature phase can be easy inves- 
tigated in the SOS model via the Wang-Landau method. 
We found that the mean-square width of the surface 
depends logarithmically on the linear system size above 
the transition temperature. This means that surface is 
rough above the faceting transition temperature. How- 
ever, this rough surface has disordered hill-and-valley 
structure and differs from typical rough phase observed 
in roughening transition (16|]. 
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